﻿//#define FaceMinCor
//#define FaceOrthCor
#define FaceOverRelax

using LxBasic;
using System;

namespace LxCore
{
    public struct FvInnFace
    {
        public FvInnFace(int[] nodes, int c0, int c1, bool isRightHand)
        {
            Nodes = new int[nodes.Length];
            Array.Copy(nodes, Nodes, nodes.Length);
            this.OSideCell = c0;
            this.NSideCell = c1;

            this.Center = new Vector3d();
            this.Normal = new Vector3d();
            this.Area = 0;
            this.OToN = new Vector3d();
            this.OToNLength = 0;
            this.gDiff = 0;
#if FaceMinCor
#else
            this.gT = 0;
#endif
            this.FaceNOrth = new();
            this.gOwner = 0;
            if (!isRightHand)
            {
                Array.Reverse(this.Nodes);
            }
        }
        public int OSideCell;
        public int NSideCell;
        public double Area;
        public double OToNLength;
        public double gOwner;//owner cell 到面的距离 投影到ownerToNeighbor，然后除以ownerToNeighbor的长度
        public double gDiff;
#if FaceMinCor
        public double gT { get => gDiff; set => gDiff = value; }//计算非正交修正时用到
#else
        public double gT;//计算非正交修正时用到
#endif
        public int[] Nodes;
        public Vector3d Center;
        public Vector3d Normal;//lenth为面的面积
        public Vector3d OToN;//由cell0 中心指向cell1中心的矢量
        public Vector3d FaceNOrth;// Normal - gT * OToN;//还要存储Normal - gT* OToN，简化计算

        public void CalGeometryInFace(Vector3d[] points)
        {
            FaceGeomCalc fg = new FaceGeomCalc();
            for (int i1 = 2; i1 < Nodes.Length; i1++)
            {
                fg.AddTriangle(points[Nodes[0]], points[Nodes[i1 - 1]], points[Nodes[i1]]);
            }
            (Normal, Center) = fg.GetResult();
            Area = Normal.Length();
        }
        public void CalGeometryAcrossInnFace(FvInnCell[] cells)
        {
            OToN = cells[NSideCell].Center - cells[OSideCell].Center;
            OToNLength = OToN.Length();
            double normalDotOToN = Vector3d.DotProduct(Normal, OToN);
#if FaceMinCor
            //此时gDiff和gT相等，gT属性访问gDiff，后面会算gT，所以这里不用算
#elif FaceOrthCor
            this.gDiff = Area / OToNLength;
#else
            this.gDiff = Normal.LengthSquared() / normalDotOToN;
#endif
            this.gT = normalDotOToN / OToN.LengthSquared();
            FaceNOrth = Normal - gT * OToN;

            Vector3d ownerToFace = Center - cells[OSideCell].Center;
            gOwner = Vector3d.DotProduct(ownerToFace, OToN) / OToN.LengthSquared();
            if (gOwner > 1) { gOwner = 1; }
            else if (gOwner < 0) { gOwner = 0; }
        }
        public double CalAverage(double[] phi)
        {
            return phi[OSideCell] * (1 - gOwner) + phi[NSideCell] * gOwner;
        }
        public Vector3d CalAverage(Vector3d[] vector)
        {
            return (1 - gOwner) * vector[OSideCell] + gOwner * vector[NSideCell];
        }
        public Vector3d CalAverage(Vector3d O_V, Vector3d N_V)
        {
            return (1 - gOwner) * O_V + gOwner * N_V;
        }
        public DiagonalMatrixR3C3 CalAverage(DiagonalMatrixR3C3[] diagMatrix)
        {
            return (1 - gOwner) * diagMatrix[OSideCell] + gOwner * diagMatrix[NSideCell];
        }
        private Vector3d calCFCorrectPart(in Vector3d aveGrad, double phiN_phiO)////使用C、F点的值对代数平均梯度修正，这是修正部分
        {
            return (phiN_phiO - Vector3d.DotProduct(aveGrad, OToN)) / (OToNLength * OToNLength) * OToN;
        }
        public Vector3d CalCFCorrectedGradientOfInner(Vector3d[] gradPhi, double[] phi)
        {
            Vector3d aveGrad = CalAverage(gradPhi);
            return aveGrad + calCFCorrectPart(aveGrad, phi[NSideCell] - phi[OSideCell]);
        }
        public Vector3d CalCFGradCorrectionOfInner(Vector3d[] gradPhi, double[] phi)//使用C、F点的值对代数平均梯度修正，获得梯度值
        {
            Vector3d aveGrad = CalAverage(gradPhi);
            return calCFCorrectPart(aveGrad, phi[NSideCell] - phi[OSideCell]);
        }

        public double CalNonOrthogonalityCorrectionOfInner(Vector3d[] gradPhi, double[] phi)
        {
#if FaceMinCor
            return Vector3d.DotProduct(CalAverage(gradPhi), FaceNOrth);
#else
            return Vector3d.DotProduct(CalAverage(gradPhi), FaceNOrth) + (phi[NSideCell] - phi[OSideCell]) * (gT - gDiff);
#endif
        }
        //public double CalNonOrthogonalityCorrection_bound(Vector3d[] gradPhi, double[] phi)
        //{
        //    return calNonOrthogonalityCorrection_FaceMinCor(gradPhi[OSideCell]);
        //}
    }
}


//public Vector3d CalONGradient_DistanceWeighted(double[] phi)//求解Grad 梯度时使用
//{
//    double deltaPhiWithDistace = (phi[NSideCell] - phi[OSideCell]) / OToNLength;
//    return deltaPhiWithDistace * OToN;
//}
//public double CalNonOrthogonalityCorrection_bound(out double a, out double b, Vector3d[] gradPhi, double[] phi,double  oCellVolumn;)//只写了匹配minCorrection的//虽然不是minCorrect也可能能用，但应该重新推导公式，进行源项线性化
//{
//    a = Vector3d.DotProduct(Normal, SmartV) / oCellVolumn;// cells[OSideCell].Volumn;//Normal - gT * OToN
//    double phiN_phiO = phi[NSideCell] - phi[OSideCell];
//    double cor = calNonOrthogonalityCorrection_FaceOrthOver(gradPhi[OSideCell], phiN_phiO);
//    b = cor - a * phiN_phiO;// (phi[NSideCell] - phi[OSideCell]);
//    return cor;
//}
//public double CalNonOrthogonalityCorrection_bound(out double a, out double b, Vector3d[] gradPhi, double[] phi, double oCellVolumn)//只写了匹配minCorrection的//虽然不是minCorrect也可能能用，但应该重新推导公式，进行源项线性化
//{
//    a = Vector3d.DotProduct(Normal, SmartV) / oCellVolumn;// cells[OSideCell].Volumn;
//    double cor = calNonOrthogonalityCorrection_FaceMinCor(gradPhi[OSideCell]);
//    b = cor - a * (phi[NSideCell] - phi[OSideCell]);
//    return cor;
//}

//        private double calNonOrthogonalityCorrection(in Vector3d meanGrad, double[] phi)
//        {
//#if FaceMinCor
//            return Vector3d.DotProduct(meanGrad, Normal - gDiff * OToN);
//#else
//            return Vector3d.DotProduct(meanGrad, Normal - _gT * OToN) + (phi[NSideCell] - phi[OSideCell]) * (_gT - gDiff);
//#endif
//        }


// a =  (Math.Pow(Area, 2) - Math.Pow(gT * OToNLength, 2)) / cells[OSideCell].Volumn;
//public double calNonOrthogonalityCorrection_Old(Vector3d[] gradPhi, double[] phi) //FVM2 P244//已废弃，用新的
//{
//    return Vector3d.DotProduct(calCFCorrectedGradient_inner(gradPhi, phi), Normal)
//            - (phi[NeighborCellIndex] - phi[OwnerCellIndex]) * gDiff;
//}
